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Abstract. A search for periodic gravitational-wave signals from isolated 
neutron stars in the NAUTILUS detector data is presented. We have analyzed 
half a year of data over the frequency band (922.2; 923.2) Hz, the spindown range 
(—1.463 X 10~*; 0) Hz/s and over the entire sky. We have divided the data into 
2 day stretches and we have analyzed each stretch coherently using matched 
filtering. We have imposed a low threshold for the optimal detection statistic 
to obtain a set of candidates that are further examined for coincidences among 
various data stretches. For some candidates we have also investigated the change 
of the signal-to-noise ratio when we increase the observation time from two to four 
days. Our analysis has not revealed any gravitational-wave signals. Therefore we 
have imposed upper limits on the dimensionless gravitational-wave amplitude over 
the parameter space that we have searched. Depending on frequency, our upper 
limit ranges from 3.4 X 10~^^ to 1.3 X 10~^^. We have attempted a statistical 
verification of the hypotheses leading to our conclusions. We estimate that our 
upper limit is accurate to within 18%. 
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1. Introduction 

We present results of the search of the NAUTILUS (a resonant bar detector pj) data 
for periodic gravitational-wave (GW) signals. We search the band (922.2; 923.2) Hz 
and the spindown range (—1.463 x 10^'*; 0)Hz/s over the entire sky. We analyze 
NAUTILUS data collected in the year 2001. We divide the data into stretches of 2 
sidereal days. Each stretch of data is analyzed coherently using matched filtering in 
the form of the J^-statistic [H [3] . We have analyzed slightly more than half a year of 
data. 

Previous analysis of bar detector data for periodic GW signals were the search of 
the galactic center and the globular cluster 47 Tucanae with the ALEGRO detector 
[4j , the search of the galactic center using the EXPLORER detector data |5j , and an 
all-sky search using the EXPLORER data [H [7] . LIGO (Laser Gravitational Wave 
Observatory) data was searched for known pulsars [H [9l [10] , over the entire sky [11] 
using the coherent method and over all the sky using incoherent methods [12j [13] . 
Currently LIGO data are analyzed over the entire sky by the Einstein@Home project 

In section[2]we present the data analysis methods used in our search. In section[3] 
we outline our search procedure. In section|4]we discuss the analysis of the candidates. 
This analysis consists of two parts: the first part is the search for coincidences between 
the candidates each obtained from a different stretch of data and the second part is an 
investigation of the increase of the signal-to-noise ratio of candidates when we increase 
the observation time from two to four days. In section [5] we impose upper limits on 
amplitudes of the gravitational waves in the parameter space that we have searched. 

2. Data analysis methods 

In order to search for gravitational waves from long lived periodic sources we have 
used the maximum likelihood (ML) method. For the case of Gaussian noise the 
ML method consists of linearly filtering the data with a template matched to the 
signal that we are searching for. The main complication of the matched filtering is 
that the signal depends on several unknown parameters. This requires evaluation 
of the likelihood function over a large parameter space. In order to minimize the 
computation time we use several data analysis tools. Firstly we find the maximum 
likelihood estimators of some parameters (4 in the case of a GW signal from a rotating 
neutron star that we are searching for) in a closed analytic form, thereby reducing the 
dimensionality of the parameter space that we have to search. The likelihood function 
over the reduced parameter space is called the J-"-statistic and it is derived in [2]. 
Secondly we analyze data of length equal to an integer multiple of a sidereal day. This 
leads to a considerable simplification of the JF-statistic and consequently reduces the 
number of numerical operations to evaluate it. The J^-statistic for observation time 
equal to an integer number of sidereal days is given in section III of }15] . Thirdly 
we use optimal numerical algorithms, in particular the Fast Fourier Transform (EFT) 
in order to calculate the J^-statistic efficiently. Fourthly we minimize the number of 
J-"-statistic calculations over the parameter space by solving a covering problem for 
this space [El [17] . Let us explain the latter two data analysis tools in more detail. 
The response of a bar detector to a gravitational-wave signal from a spinning neutron 
star is summarized in section 2.1 of [7|. 



All-sky search of NAUTILUS data 



3 



Fast Fourier Transform. Estimates have shown [211 120] that for the bandwidth 
and the spindown range that we search we need to take into account in our templates 
only one spindown parameter in order to match the signal. Consequently the phase 
modulation function (f>(t) of the waveform is given by 

^ ujot + Loit^ + {loo + 2cjit) "° ' '^'^'-^^ , (1) 

c 

where ujq is angular frequency and uji is the spindown parameter, ng is the constant 
unit vector in the direction of the star in the Solar System Barycenter (SSB) reference 
frame (it depends on the right ascension a and the declination 6 of the source) , and 
Td is the vector joining the SSB with the detector and c is the speed of light. The 
detection statistic T involves two integrals of the form 

F= r" x{t) a{t) e-'*^'Ut, (2) 
Jo 

where x{t) is the data stream, a{t) is the amplitude modulation function that depends 
on 5 and a. The above integral is not a Fourier transform because the frequency 
coq in the phase multiplies the term no • rd(^) which is a non linear function of time. 
In order to convert the integral into a Fourier transform we introduce the following 
interpolation procedure. The phase (f){t) [Equation ([1])] can be written as 

(/)(t) =tJo[t + 0mW] (3) 



where 



M---^^^^, (4) 



^2 , o^O-I-dW^ 



(jjsit) ■.= ujit^ + 2— — ^^ujit. (5) 
c 

The functions 0m (0 0s (0 depend on the angular frequency loq. We can 

write the integral ^ as 

[ ° x{t) a{t) e-^'>=(*) exp { - iujo[t + 0m (t)]} dt. (6) 



We see that the integral ([6]) can be interpreted as a Fourier transform (and computed 
efficiently with an FFT), if 0m = 0. In order to convert equation ([5]) to a Fourier 
transform we introduce a new time variable tj,, so called barycentric time \TE[ f^. 

4:=t + 0m(t). (7) 

In the new time coordinate the integral ^ is approximately given by (see Ref. [2], 
Sec. HID) 

F^ /^\-[t(4)]a[i(4)]e-*^=[*(*^)le-''"«*''d4. (8) 



Thus in order to compute the integral ([6]), we first multiply the data x{t) by the 
function a(t) exp[— i0s(i)] for each set of the parameters u!i,6,a and then resample 
the resulting function according to equation d?]). At the end we perform the FFT. 

The covering problem. The covering problem is to find the minimum number of 
templates in the parameter space |17| . so that the fractional loss in signal to ratio is 
not less than 1 — MM {MM is the minimal match parameter introduced by Owen 
|22]). In order to solve the covering problem we introduce a useful approximate 
model of the gravitational-wave signal from a rotating neutron star. The model 
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relies on (i) neglecting all spindowns in the phase modulation due to motion of the 
detector with respect to the SSB; and (ii) discarding the component of the vector 
(connecting the SSB and the detector) which is perpendicular to the ecliptic plane. 
These approximations lead to the following phase model of the signal: 

(/)ii„(i) = LOot + uit^ + aifii{t) + a2M2(0: (9) 

where ai and a2 are new constant parameters, 

ai Wo (sin a cos (5 cose + sin (5 sine), (10) 

a2 :— ujQCOsacoaS, (11) 

where e is the obliquity of the ecliptic and where /ii(t) and fi2(t) are known functions 
of time, 

:=^EsW + ^EWcose, (12) 

^i2{t):=RUt) + Rm- (13) 
i?gg is the x-component of the vector joining the center of Earth and the SSB, and i?| 
is the cc-component of the vector joining the center of Earth and the detector. R^^it) 
and R^{t) are the corresponding y-components. We also neglect the slowly varying 
modulation of the signal's amplitude, so finally we approximate the whole signal h(t) 

by 

h{t)^ Ao cos {(j)iinit) + (j)o), (14) 

where Aq and (po are the constant amplitude and initial phase, respectively. The 
above signal model is called linear because it has the property that its phase given 
by equation ([5]) is a linear function of the parameters. We have shown [TH] that the 
above model is a good approximation to the accurate response of the detector to the 
GW signal in the sense that the Fisher matrix for the linear model reproduces well the 
Fisher matrix for the accurate model. Thus whenever a Fisher matrix is needed we can 
use the Fisher matrix for the linear model as an approximation to the Fisher matrix 
for the accurate model. The great advantage of the linear model is that components 
of its Fisher matrix are constant, independent of the values of the parameters. In 
order to solve the covering problem for the parameter space we use the Fisher matrix 
as a metric on the parameter space. Because the components of the Fisher matrix are 
constant the grid is uniform what greatly simplifies its construction. In our search, as 
a grid we use the hypercubic lattice [16j . However we have an additional constraint. 
In order to apply the EFT algorithm the nodes of the grid must coincide with the 
Fourier frequencies. We have constructed a suitable grid by performing rotations and 
dilatations of the original hypercubic lattice. We construct the grid in the parameters 
LOojUJi^ai, oi2 and then transform it to parameters wo, ^i, for which the ^-statistic 
is calculated. 

The linear parametrization has one more application. We use it in order to 
calculate the threshold for the ^-statistic corresponding to a certain false alarm 
probability. Namely, using the linear parametrization we divide the parameter space 
into cells as explained in |2l [20] . All the cells are exactly the same and their number 
A^c is easily calculated using the Fisher matrix (see section IIIB of [20]). The false 
alarm probability a is the probability that T exceeds threshold To in one or more 
cells and is given by 

a^\-{\-PF{To)X'\ (15) 
where Pp is the false alarm probability for a single cell. 
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3. Search procedure 

We have searched the data collected by the NAUTILUS detector in the year 2001. 
The bandwidth of (922.2; 923.2) Hz, where the detector is most sensitive, has been 
analyzed. We have divided the data into stretches which span two sidereal days. 
We have assumed a minimum pulsar spindown age Tmin equal to 1000 yrs and so 
we have searched the negative frequency time derivatives in the range of (—1.463 x 
10~^; 0) Hz/s. For this Tmin and two days of the observation time it is sufhcient to 
include only one spindown in the phase [HI |^ . Each two day sequence was analyzed 
coherently using the J^-statistic. We have used the constrained hypercubic grid as 
explained in the previous section. For the grid construction we have assumed the 
minimal match parameter MM — V3/2 [22j . Using this minimal-match value our grid 
consists of around 3.1 x 10^^ grid points (2^^ frequency bins, ~ 10^ spindowns, ~ 6x 10"* 
sky positions). The threshold on 2!F corresponding to 1% false alarm probability has 
been calculated using Equation (jisp and is around 72. In order to compensate the loss 
of signal-to-noise ratio (SNR) due to the discreteness of the grid, imperfect templates 
and numerical approximation in evaluation of the jT-statistic (resampling procedure) 
we have adopted two lower thresholds on 2!F equal to 40 and 50. We have registered 
parameters of all templates which crossed the threshold of 40. For threshold crossings 
of 50 we have performed a verification procedure. The verification procedure consisted 
of calculating the JF-statistic for the template parameters of the candidate using a four 
day stretch of data involving the original two day stretch. For a true gravitational- 
wave signal by this procedure one would expect an increase of signal-to-noise ratio 
by y/2. In total, we have analyzed 93 two day data stretches. In Figure [1] we have 
presented the two-sided amplitude spectrum of the NAUTILUS detector data that we 
have analyzed. The spectrum was obtained in the following way. We have estimated 
the power spectrum density in each of the 93 two day data sequences and then we 
have taken the square root of the average of the 93 power spectra. We see that the 
best sensitivity is around 5 x 10~^^ Hz~^/^. Moreover we have obtained the rms error 
of our power spectrum estimate by calculating the variance from the estimates of the 
spectra of in each of the 93 data segments. The relative 1 a error in the amplitude 
power spectrum is around 18%. 

During the search we have obtained 537 665 380 candidates above the 2J^- 
threshold of 40 and 9 038 817 above the 2J^-threshold of 50. 

4. Analysis of the candidates 

4.1- Signal-to-noise ratio of the candidates 

In Figure [2] we have plotted a histogram of the frequencies of all the candidates 
above the 2J^-threshold of 50. The histogram shows an excess of candidates in the 
frequency band of (922.4; 922.6) Hz. This excess is a result of the presence of a periodic 
interference in the data that appears as a series of harmonics in the bandwidth of the 
detector. One of the harmonics is located in the subband (922.2; 923.2) Hz. The effect 
of the harmonic is visible in our estimate of the spectrum (Figure [T|) and appears as 
a bump in the band (922.5; 922.6) Hz. 

As a first step in the candidate analysis we have calculated the increase in signal- 
to-noise ratio when we increase the observation time from two to four days. This 
has been done for all the candidates above the threshold ~ 50. Figure [3] shows 




922.2 922.4 922.6 922.8 923 923.2 

Frequency [Hz] 

Figure 1. Estimation of tiic two-sided amplitude spectrum of the NAUTILUS 
data in the year 2001 and the rms error of the estimate. The thick hne shows the 
estimate and the two thin hnes correspond to the 1 a error. 



the highest increase in SNR for candidates when going from a two day data stretch 
to the four day one. The maximum is calculated for each of the 93 data stretches 
analyzed. We see that typically the highest gain in the signal-to-noise is 1.2. This 
should be compared with the theoretical gain of V2 of SNR when we increase the 
observation time by a factor of 2. The periodic interference present in the data to 
which we attribute these maximum SNR increases does not gives a higher increase of 
the SNR because its frequency changes erratically over the observation time of days 
and it cannot reproduce the Doppler shift of a real GW signal modulated by detector 
motion with respect to the SSB. Assuming that the two day sequence is independent 
of the four day sequence we could perform the F-test that consists of calculation of 
the ratio F of the jF-statistic for 4 days observation time and the jF-statistic for 2 days 
observation time. Taking as the null hypothesis for the test that data is only Gaussian 
noise the 2^-statistic has the central distribution with 4-degrees of freedom and the 
ratio F has Fisher-Snedeckor distribution F(4, 4). The typical highest value of F for 
a given data segment is around 1.5 The probability of F crossing the threshold 1.5 is 
around 37%. This would give a high confidence that data is noise only. Unfortunately 
this is only a crude approximation because the two day sequence is contained in the 
four day one and the assumption of independence of the two JF-statistic is not fulfilled. 
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Figure 2. Histogram of the frequencies of candidates obtained in the search of 
all 93 two day data stretches above a 2^-threshoId of 50. 



4-2. Coincidences among the candidates from different data stretches 

Candidates from different data stretches are considered coincident if they cluster 
closely together in the four-dimensional parameter space {ujQ,uji,S,a). We employ 
the clustering method described in [13], which uses a grid of "coincidence cells". 
This method will reliably detect strong signals which would produce candidates with 
closely-matched parameters in many of the different data stretches. 

In a first step, the frequency value of each candidate above the threshold of 
2!F = 40 is shifted to the same fiducial time: the GPS start time of the earliest (j — 1) 
stretch, ^fiducial = = 662 547 735.9988098 s. Defining Tq to be the time span of two 
sidereal days, the frequencies of the candidates are shifted to ^fiducial via 

t^o(ifiduciai) = ^o{tj) + {j - 1) 2cJiTo, (16) 

where tj is the starting time of the j'th data stretch, given by tj — ^fiducial + (j ^ 1) ^q. 

To find coincidences, a grid of cells is constructed such that the cells are 
rectangular in the coordinates {u>o,uJi,S,acos6). The dimensions of the cells are 
adapted to the parameter space search. Thus, the cells are constructed to be as small 
as possible to reduce the probability of coincidences due to false alarms. However, 
since each of the 93 different data stretches uses a slightly different parameter space 
grid, the coincidence cells must be chosen to be large enough that the candidates from 
a source (which would appear at slightly different points in parameter space in each 
of the 93 data stretches) would still lie in the same coincidence cell. As a conservative 
choice we use cell sizes in uiq of 5.8 x lO^^'Hz, in ui of 2.08 x 10^^"'^ Hz s^"'^, and an 
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Figure 3. Highest increase (vertical axis) in signal-to-noise ratio for candidates 
in each of the 93 data stretches analyzed. The two day stretches of Nautilus 2001 
data are numbered form 1 to 182. The missing lines in the plot indicate that the 
corresponding data stretch was not analyzed. 



isotropic cell grid in the sky with equatorial spacing of 0.028 rad. Each candidate 
event is assigned to a particular cell. In cases where two or more candidate events 
from the same data stretch j fall into the same cell, only the candidate having the 
largest value of 2J^ is retained in the cell. Then the number of candidate events per 
cell coming from distinct data stretches is counted. 

From the 93 different data stretches, this coincidence method found that we get 
candidates which appear consistently in no more than 4 data stretches uniformly over 
the search bandwidth, where there are no instrumental interferences. This is the 
background of the number of coincidences. We would like to test the null hypothesis 
that the coincidences are result of the noise only. Let us assume that the parameter 
space is divided into Nceii independent coincidence cells, the candidate events are 
independent and the probability for a candidate event to fall into any given coincidence 
cell is 1 = 1 /Nceii ■ Thus probability e that a given coincidence cell is populated with 
one or more candidate event is given by 

6 = i-(i--^r-, (17) 

where Sseg is the number of candidate events per data segment. The probability pp 
that any given coincidence cell contains candidate events from €„ 
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data segments is given by a binomial distribution 

P^= E (18) 

Finally the probability Pp that there is Cmax or more coincidences in one or more of 
the Nceii cells is 

Pf = 1 - {1 ~ Pf)''"" ■ (19) 

The expected number of cells with Cmax or more coincidences is given by 

Nf^N.mpf. (20) 

In our case the number of cells is given by N^eii = 5-9 X lO^", the number of 
data segments is Ngeg = 93, and the number of candidates per data segment is 
£seg = 5.8 X 10^. From Equation we find that the probability of finding Cmax = 4 
or more coincident candidates is almost one. Thus for the background coincidences we 
can accept the null hypothesis that they are from noise only with a high confidence. 
Over the bandwidth (922.4; 922.6) Hz we find an excess of coincidences with the 
maximum of 8 coincidences. By Equation p^ . the false alarm probability associated 
with with 8 or more coincidences is of the order of 10~^^ and thus they cannot 
be attributed to noise. We consider these coincidences to be due to the periodic 
interference present in the data. 



5. Upper limits 

Our verification procedure consisting of coincidences among the candidates from 
distinct data segments and an analysis of the increase of signal-to-noise ratio presented 
in section |4] did not produce convincing evidence of a gravitational-wave signal. 
We therefore proceeded to estimate the upper limits for the amplitudes of the 
gravitational- wave signals in the parameter space that we have searched. Detection of 
a signal is signified by a large value of the jT-statistic that is unlikely to arise from the 
noise-only distribution. If instead the value of is consistent with pure noise with 
high probability we can place an upper limit on the strength of the signal. One way 
of doing this is to take the loudest event obtained in the search and solve the equation 

P =PD(pul,.Floudcst) (21) 

for signal-to- noise ratio pui, where Pd is the detection probability, J^oudest is the value 
of the ^-statistic corresponding to the loudest event, and P is a chosen confidence. 
Then is the desired upper limit with confidence V. We can also obtain an upper 
limit with confidence V for several independent searches from the equation 

L 

V = l-Y[[l- PoiPuh^Floudcsts)], (22) 
s=l 

where Jioudcsts is the threshold corresponding to the loudest event in s'th search and 
L is the number of searches. Here V is the probability that a signal with signal-to- 
noise ratio Pui crosses the threshold .Fioudosts in at least one of the L independent 
searches. To calculate p^i we assume that the data have a Gaussian distribution and 
consequently the probability of detection has a non-central distribution with 4 
degrees of freedom and the noncentrality parameter equal to . We have investigated 
this assumption by obtaining histograms of the 2J^-statistic values of the candidates 




Figure 4. Probability distribution of 2.7-"-statistic values of the candidates. The 
light lines are obtained from histograms of the 2T values of each data segment. 
The thick line represents the theoretical central distribution with 4 degrees of 
freedom. 



and comparing them to the central distribution with 4 degrees of freedom. The 
result is shown in Figure There is an overall qualitative agreement of candidates 
distributions with the theoretical one. However, the candidates distributions do not 
pass a goodness-of-fit test for a distribution at the significance level of 5%. 

In order to translate our upper limit on the SNR into the upper limit on the 
gravitational-wave amplitude, we use the Equation (93) of [2] for signal-to-noise ratio 
of a GW signal from a spinning neutron star averaged over the source position and 
orientation. Thus /lui and pui are related by the following formula: 

hM)^l^^Pni, (23) 

where S{f) is one-sided spectral density at frequency /. We have used Equations ([22| 
and (|23p to obtain upper limits in 0.1 Hz bands over the bandwidth (922.2; 923.2) Hz 
that we have searched. The upper limit results are presented in Figure [S] Assumed 
Gaussian noise, we have chosen the confidence V = 90% and we denote the upper 
limits by h^^'^". Our best upper limit is equal to 3.4 x 10^^'^ at a frequency of 922.55 
Hz. Using our 1 a rms error of the amplitude power spectrum estimate we reckon that 
our upper limit has likewise an error of 18%. 
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Figure 5. Upper limits based on the loudest candidate for 0.1 Hz frequency 
bands over the 1 Hz bandwidth searched. 
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